//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2012 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : boundarypsi_kernels.h
// Author     : M.A. de Jong
//----------------------------------------------------------------------------

#ifndef BOUNDARYPSI_KERNELS_H
#define BOUNDARYPSI_KERNELS_H

namespace wavescomputer
{

namespace boundarypsi
{

namespace reflectiveneumann
{

// Reflective Neumann boundary (North version)
template <class T>
__global__ void kernel_bp_rn_north(T* psi, 
                                   const T* psiN,
                                   const T* psi_in_N_ext, 
                                   const T* psi_in_N_cen, 
                                   const int ld,
                                   const int Nx1,
                                   const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH + tid) * ld + (BORDER_WIDTH + Nx2); // not coalesced

    if (tid < Nx1) 
    {
        psi[tid2D] = psi[tid2D - 1] + psi_in_N_ext[tid1D] - psi_in_N_cen[tid1D];
    }
}

// Reflective Neumann boundary (East version)
template <class T>
__global__ void kernel_bp_rn_east(T* psi, 
                                  const T* psiE,
                                  const T* psi_in_E_ext, 
                                  const T* psi_in_E_cen, 
                                  const int ld,
                                  const int Nx1,
                                  const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH + Nx1) * ld + BORDER_WIDTH + tid;   // coalesced

    if (tid < Nx2) 
    {
        psi[tid2D] = psi[tid2D - ld] + psi_in_E_ext[tid1D] - psi_in_E_cen[tid1D];
    }
}

// Reflective Neumann boundary (South version)
template <class T>
__global__ void kernel_bp_rn_south(T* psi, 
                                   const T* psiS,
                                   const T* psi_in_S_ext, 
                                   const T* psi_in_S_cen, 
                                   const int ld,
                                   const int Nx1,
                                   const int Nx2) 
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH + tid) * ld + BORDER_WIDTH - 1;     // not coalesced

    if (tid < Nx1) 
    {
        psi[tid2D] = psi[tid2D + 1] + psi_in_S_ext[tid1D] - psi_in_S_cen[tid1D];
    }
}

// Reflective Neumann boundary (West version)
template <class T>
__global__ void kernel_bp_rn_west(T* psi, 
                                  const T* psiW,
                                  const T* psi_in_W_ext, 
                                  const T* psi_in_W_cen, 
                                  const int ld,
                                  const int Nx1, 
                                  const int Nx2) 
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH - 1) * ld + BORDER_WIDTH + tid;     // coalesced

    if (tid < Nx2) 
    {
        psi[tid2D] = psi[tid2D + ld] + psi_in_W_ext[tid1D] - psi_in_W_cen[tid1D];
    }
}

} // namespace reflectiveneumann

namespace weaklyreflective
{

// weakly-reflective boundary (North version)
template <class T>
__global__ void kernel_bp_wr_north(T* psi, 
                                   const T* psiN, 
                                   const T* psi_in_N_ext, 
                                   const int ld,
                                   const int Nx1,
                                   const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH + tid) * ld + (BORDER_WIDTH + Nx2); // not coalesced

    if (tid < Nx1) 
    {
        psi[tid2D] = psiN[tid1D] + psi_in_N_ext[tid1D];
    }
}

// weakly-reflective boundary (East version)
template <class T>
__global__ void kernel_bp_wr_east(T* psi, 
                                  const T* psiW, 
                                  const T* psi_in_W_ext, 
                                  const int ld,
                                  const int Nx1,
                                  const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH + Nx1) * ld + BORDER_WIDTH + tid;   // coalesced

    if (tid < Nx2) 
    {
        psi[tid2D] = psiW[tid1D] + psi_in_W_ext[tid1D];
    }
}

// weakly-reflective boundary (South version)
template <class T>
__global__ void kernel_bp_wr_south(T* psi, 
                                   const T* psiS, 
                                   const T* psi_in_S_ext, 
                                   const int ld,
                                   const int Nx1,
                                   const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH + tid) * ld + BORDER_WIDTH - 1;     // not coalesced

    if (tid < Nx1) 
    {
        psi[tid2D] = psiS[tid1D] + psi_in_S_ext[tid1D];
    }
}

// weakly-reflective boundary (West version)
template <class T>
__global__ void kernel_bp_wr_west(T* psi, 
                                  const T* psiW, 
                                  const T* psi_in_W_ext, 
                                  const int ld,          // leading dimension
                                  const int Nx1,         // no. nodes in y-direction
                                  const int Nx2)         // no. nodes in x-direction
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH - 1) * ld + BORDER_WIDTH + tid;     // coalesced

    if (tid < Nx2) 
    {
        psi[tid2D] = psiW[tid1D] + psi_in_W_ext[tid1D];
    }
}

} // namespace weaklyreflective

} // namespace boundarypsi

} // namespace wavescomputer

#endif // BOUNDARYPSI_KERNELS_H

